###############################################################################################
# Mean Comparisons
data$inter <- interaction(data$rel_country, data$rel_devotion)
data$inter <- interaction(data$inter, data$rel_religion)
data$inter <- interaction(data$inter, data$rel_banperm)

labels <- c("Bulgaria.devout.Christians.ban", 
            "Nigeria.devout.Christians.ban", 
            "native.devout.Christians.ban", 
            "Bulgaria.non-practising.Christians.ban", 
            "Nigeria.non-practising.Christians.ban", 
            "native.non-practising.Christians.ban", 
            "Bulgaria.radical.Christians.ban", 
            "Nigeria.radical.Christians.ban", 
            "native.radical.Christians.ban", 
            "Bulgaria.devout.Muslims.ban", 
            "Nigeria.devout.Muslims.ban", 
            "native.devout.Muslims.ban", 
            "Bulgaria.non-practising.Muslims.ban", 
            "Nigeria.non-practising.Muslims.ban",
            "native.non-practising.Muslims.ban", 
            "Bulgaria.radical.Muslims.ban",
            "Nigeria.radical.Muslims.ban", 
            "native.radical.Muslims.ban", 
            "Bulgaria.devout.Christians.permit", 
            "Nigeria.devout.Christians.permit",  
            "native.devout.Christians.permit", 
            "Bulgaria.non-practising.Christians.permit", 
            "Nigeria.non-practising.Christians.permit", 
            "native.non-practising.Christians.permit", 
            "Bulgaria.radical.Christians.permit", 
            "Nigeria.radical.Christians.permit", 
            "native.radical.Christians.permit", 
            "Bulgaria.devout.Muslims.permit", 
            "Nigeria.devout.Muslims.permit", 
            "native.devout.Muslims.permit",
            "Bulgaria.non-practising.Muslims.permit", 
            "Nigeria.non-practising.Muslims.permit", 
            "native.non-practising.Muslims.permit", 
            "Bulgaria.radical.Muslims.permit", 
            "Nigeria.radical.Muslims.permit", 
            "native.radical.Muslims.permit")


##########################################################################################
# Feeling Score

feel.means <- aggregate(data$rel_feel, by=list(g=data$inter), mean, na.rm=T)
feel.means[,2] <- round(feel.means[,2])
rownames(feel.means) <- labels

# natives
nat.non <- feel.means[seq(3, 36, 3),][seq(2, 12, 3),]
nat.dev <- feel.means[seq(3, 36, 3),][seq(1, 12, 3),]
nat.rad <- feel.means[seq(3, 36, 3),][seq(3, 12, 3),]

nat.non <- cbind(nat.non[1:2,2], nat.non[3:4,2])
nat.dev <- cbind(nat.dev[1:2,2], nat.dev[3:4,2])
nat.rad <- cbind(nat.rad[1:2,2], nat.rad[3:4,2])

# bulgaria
bul.non <- feel.means[seq(1, 36, 3),][seq(2, 12, 3),]
bul.dev <- feel.means[seq(1, 36, 3),][seq(1, 12, 3),]
bul.rad <- feel.means[seq(1, 36, 3),][seq(3, 12, 3),]

bul.non <- cbind(bul.non[1:2,2], bul.non[3:4,2])
bul.dev <- cbind(bul.dev[1:2,2], bul.dev[3:4,2])
bul.rad <- cbind(bul.rad[1:2,2], bul.rad[3:4,2])

# nigeria
nig.non <- feel.means[seq(2, 36, 3),][seq(2, 12, 3),]
nig.dev <- feel.means[seq(2, 36, 3),][seq(1, 12, 3),]
nig.rad <- feel.means[seq(2, 36, 3),][seq(3, 12, 3),]

nig.non <- cbind(nig.non[1:2,2], nig.non[3:4,2])
nig.dev <- cbind(nig.dev[1:2,2], nig.dev[3:4,2])
nig.rad <- cbind(nig.rad[1:2,2], nig.rad[3:4,2])

non <- rbind(nat.non, nat.non[2,]-nat.non[1,], bul.non, bul.non[2,]-bul.non[1,], nig.non, nig.non[2,]-nig.non[1,])
dev <- rbind(nat.dev, nat.dev[2,]-nat.dev[1,], bul.dev, bul.dev[2,]-bul.dev[1,], nig.dev, nig.dev[2,]-nig.dev[1,])
rad <- rbind(nat.rad, nat.rad[2,]-nat.rad[1,], bul.rad, bul.rad[2,]-bul.rad[1,], nig.rad, nig.rad[2,]-nig.rad[1,])

non.diff <- non[,2]-non[,1]
dev.diff <- dev[,2]-dev[,1]
rad.diff <- rad[,2]-rad[,1]
non.diff[c(3,6,9)] <- NA
dev.diff[c(3,6,9)] <- NA
rad.diff[c(3,6,9)] <- NA

# Mean Differences
# across policy
p.s <- rep(NA, 18)
for(i in 1:18){
  p.s[i] <- t.test(data$rel_feel[data$inter==feel.means[i,1]], data$rel_feel[data$inter==feel.means[i+18,1]])$p.value 
}

p.non <- p.s[c(4:6, 13:15)]*45
p.dev <- p.s[c(1:3, 10:12)]*45
p.rad <- p.s[c(7:9, 16:18)]*45

p.non <- round(c(p.non[3], p.non[6], NA, p.non[1], p.non[4], NA, p.non[2], p.non[5], NA), 2)
p.dev <- round(c(p.dev[3], p.dev[6], NA, p.dev[1], p.dev[4], NA, p.dev[2], p.dev[5], NA), 2)
p.rad <- round(c(p.rad[3], p.rad[6], NA, p.rad[1], p.rad[4], NA, p.rad[2], p.rad[5], NA), 2)

feel.tab <- cbind(non, non.diff, p.non, dev, dev.diff, p.dev, rad, rad.diff, p.rad)
colnames(feel.tab) <- rep(c("Ban", "Permit", "Diff", "p"), 3)

# across religion
# ban
p.s.a <- rep(NA, 9)
for(i in 1:9){
  p.s.a[i] <- t.test(data$rel_feel[data$inter==feel.means[i,1]], data$rel_feel[data$inter==feel.means[i+9,1]])$p.value 
}

# perm
p.s.b <- rep(NA, 9)
for(i in 19:27){
  p.s.b[i-18] <- t.test(data$rel_feel[data$inter==feel.means[i,1]], data$rel_feel[data$inter==feel.means[i+9,1]])$p.value 
}

p.nat <- cbind(p.s.a, p.s.b)[seq(3, 9, 3),]
p.bul <- cbind(p.s.a, p.s.b)[seq(1, 9, 3),]
p.nig <- cbind(p.s.a, p.s.b)[seq(2, 9, 3),]

p.nat <- c(p.nat[2,], NA, NA, p.nat[1,], NA, NA, p.nat[3,], NA, NA)
p.bul <- c(p.bul[2,], NA, NA, p.bul[1,], NA, NA, p.bul[3,], NA, NA)
p.nig <- c(p.nig[2,], NA, NA, p.nig[1,], NA, NA, p.nig[3,], NA, NA)

p.nat <- round(p.nat*45, 2)
p.bul <- round(p.bul*45, 2)
p.nig <- round(p.nig*45, 2)

feel.tab <- rbind(feel.tab[1:3,], p.nat, feel.tab[4:6,], p.bul, feel.tab[7:9,], p.nig)
rownames(feel.tab) <- rep(c("Christian", "Muslim", "Diff", "p"), 3)


################################################################################################
# Support
data$support <- NA
data$support <- ifelse(data$rel_banperm==1, as.numeric(data$rel_posban), NA)
data$support <- ifelse(data$rel_banperm==2, as.numeric(data$rel_posper), data$support)
data$support <- data$support * -1 + 6

support.means <- aggregate(data$support, by=list(g=data$inter), mean, na.rm=T)
support.means[,2] <- round(support.means[,2],1)
rownames(support.means) <- labels

# natives
nat.non <- support.means[seq(3, 36, 3),][seq(2, 12, 3),]
nat.dev <- support.means[seq(3, 36, 3),][seq(1, 12, 3),]
nat.rad <- support.means[seq(3, 36, 3),][seq(3, 12, 3),]

nat.non <- cbind(nat.non[1:2,2], nat.non[3:4,2])
nat.dev <- cbind(nat.dev[1:2,2], nat.dev[3:4,2])
nat.rad <- cbind(nat.rad[1:2,2], nat.rad[3:4,2])

# bulgaria
bul.non <- support.means[seq(1, 36, 3),][seq(2, 12, 3),]
bul.dev <- support.means[seq(1, 36, 3),][seq(1, 12, 3),]
bul.rad <- support.means[seq(1, 36, 3),][seq(3, 12, 3),]

bul.non <- cbind(bul.non[1:2,2], bul.non[3:4,2])
bul.dev <- cbind(bul.dev[1:2,2], bul.dev[3:4,2])
bul.rad <- cbind(bul.rad[1:2,2], bul.rad[3:4,2])

# nigeria
nig.non <- support.means[seq(2, 36, 3),][seq(2, 12, 3),]
nig.dev <- support.means[seq(2, 36, 3),][seq(1, 12, 3),]
nig.rad <- support.means[seq(2, 36, 3),][seq(3, 12, 3),]

nig.non <- cbind(nig.non[1:2,2], nig.non[3:4,2])
nig.dev <- cbind(nig.dev[1:2,2], nig.dev[3:4,2])
nig.rad <- cbind(nig.rad[1:2,2], nig.rad[3:4,2])

non <- rbind(nat.non, nat.non[2,]-nat.non[1,], bul.non, bul.non[2,]-bul.non[1,], nig.non, nig.non[2,]-nig.non[1,])
dev <- rbind(nat.dev, nat.dev[2,]-nat.dev[1,], bul.dev, bul.dev[2,]-bul.dev[1,], nig.dev, nig.dev[2,]-nig.dev[1,])
rad <- rbind(nat.rad, nat.rad[2,]-nat.rad[1,], bul.rad, bul.rad[2,]-bul.rad[1,], nig.rad, nig.rad[2,]-nig.rad[1,])

non.diff <- non[,2]-non[,1]
dev.diff <- dev[,2]-dev[,1]
rad.diff <- rad[,2]-rad[,1]
non.diff[c(3,6,9)] <- NA
dev.diff[c(3,6,9)] <- NA
rad.diff[c(3,6,9)] <- NA

# Mean Differences
# across policy
# support
p.s <- rep(NA, 18)
for(i in 1:18){
  p.s[i] <- t.test(data$support[data$inter==support.means[i,1]], data$support[data$inter==support.means[i+18,1]])$p.value 
}

p.non <- p.s[c(4:6, 13:15)]*45
p.dev <- p.s[c(1:3, 10:12)]*45
p.rad <- p.s[c(7:9, 16:18)]*45

p.non <- round(c(p.non[3], p.non[6], NA, p.non[1], p.non[4], NA, p.non[2], p.non[5], NA), 2)
p.dev <- round(c(p.dev[3], p.dev[6], NA, p.dev[1], p.dev[4], NA, p.dev[2], p.dev[5], NA), 2)
p.rad <- round(c(p.rad[3], p.rad[6], NA, p.rad[1], p.rad[4], NA, p.rad[2], p.rad[5], NA), 2)

support.tab <- cbind(non, non.diff, p.non, dev, dev.diff, p.dev, rad, rad.diff, p.rad)
colnames(support.tab) <- rep(c("Ban", "Permit", "Diff", "p"), 3)

# across religion
# ban
p.s.a <- rep(NA, 9)
for(i in 1:9){
  p.s.a[i] <- t.test(data$support[data$inter==support.means[i,1]], data$support[data$inter==support.means[i+9,1]])$p.value 
}

# perm
p.s.b <- rep(NA, 9)
for(i in 19:27){
  p.s.b[i-18] <- t.test(data$support[data$inter==support.means[i,1]], data$support[data$inter==support.means[i+9,1]])$p.value 
}

p.nat <- cbind(p.s.a, p.s.b)[seq(3, 9, 3),]
p.bul <- cbind(p.s.a, p.s.b)[seq(1, 9, 3),]
p.nig <- cbind(p.s.a, p.s.b)[seq(2, 9, 3),]

p.nat <- c(p.nat[2,], NA, NA, p.nat[1,], NA, NA, p.nat[3,], NA, NA)
p.bul <- c(p.bul[2,], NA, NA, p.bul[1,], NA, NA, p.bul[3,], NA, NA)
p.nig <- c(p.nig[2,], NA, NA, p.nig[1,], NA, NA, p.nig[3,], NA, NA)

p.nat <- round(p.nat*45, 2)
p.bul <- round(p.bul*45, 2)
p.nig <- round(p.nig*45, 2)

support.tab <- rbind(support.tab[1:3,], p.nat, support.tab[4:6,], p.bul, support.tab[7:9,], p.nig)
rownames(support.tab) <- rep(c("Christian", "Muslim", "Diff", "p"), 3)

################################################################################################
